The X-Philes (2nd Revision)/The X-Philes Number 1 (1995).iso/xphiles/hp48_1/polyfit1

Internet Message Format | 1995-03-31 | 12KB

Path: seq!spell
From: David F. Kurth <[e-mail]>
Subject:  v01i011:  polyfit10 - Polyfit v1.0: produce polynomial given N data points, Part01/01
Newsgroups: comp.sources.hp48
Followup-To: comp.sys.hp48
Approved: [e-mail]

Checksum: 3154341908 (verify with brik -cv)

Submitted-by: David F. Kurth <[e-mail]>
Posting-number: Volume 1, Issue 11
Archive-name: polyfit10/part01

BEGIN_DOC polyfit.doc
I've always thought there should be a good use for a program like this,
but I still haven't found it.  It was a good challenge to program the
algorithm, and maybe someone else can find a good use for it.

POLYFIT produces the coefficients of an N-1 degree polynomial given
N data points that lie on the curve of that polynomial.
For example, consider the following table of data points:

   X     Y
-----------
  -2     1
  -1    -1
   0     1
   1    -1
   2     1

There are 5 data points, and the Y value a very symmetric.  POLYFIT
produces the following polynomial:

  f(x) = 2/3*x^4 + 0*x^3 - 8/3*x^2 + 0*x^1 + 1*x^0

If you plot this with an X range between -2 and 2, you will see the
curve "oscillate" between +1 and -1 as it hits all five data points.

POLYFIT is not a least square's type curve fit program.  There is no
minimization of errors.  The polynomial coefficients found by POLYFIT
will exactly hit the given data points (to within the accuracy of the
calculator).

Using POLYFIT
Enter the following code into your calculator.  After loading, enter
the new directory if not already there.  Push the CST key load the
custom menu which should reveal the following keys:

  POLY - Executes the main program which takes the input array
         of X-Y points from the stack and produces the K array of
         polymonial coefficients.
  VEIW - Displays the coefficients of the resulting polynomial.
  FX   - Calculates f(x) using level 1 of the stack for x.
  X-Y  - Recalls the input data points to the stack for editing
         or viewing.

To begin, you must first enter the input data points array.  This
is easiest to do by using the Matrix Writer (Right-shift ENTER).
Enter your first X-Y pair of points, then push the down arrow
which informs the calculator that the array is only 2 elements wide.
Continue entering the rest of your points until done.  Then just
push ENTER to return the array to the stack.  Now push POLY to
calculate the coefficients.

The display should produce "DOING L:" which can be ignored (the
program calculates a triangular array of difference values used
for the rest of the computations).  After this, the Ki coefficient
values are calculated and briefly displayed.

When done, the K values can be displayed again by using the VIEW
key.  Output values can be calculated by entering an X value on
the stack and pressing the FX key.  Plotting the equation can
be done by pressing Left-shift PLOT, PLOTR, and then AUTO.

For a small number of data points, the calculations are fairly
quick.  For each additional data point entered, calculation time
roughly doubles.  I've calculated 15 data points, but it took the
calculator about 7.5 hours to finish.  However, with sufficient
memory and batteries, have at it.

Dave Kurth
[e-mail]

END_DOC

BYTES: #B610h 2267.5

BEGIN_RPL polyfit
%%HP: T(3)A(D)F(.);
DIR
POLY
@ This program takes an array of N X-Y values from the stack @
@ and calculates the coefficients of an N-1 polynomial that @
@ intersect with all N input points.  The N X-Y values are most @
@ easily entered using the Matrix Writer Application. @
@@
@ Directory @
@ Checksum # B610h
@ Bytes    2267.5
@@
\<< IF DEPTH                            @ Check for an object on stack
  THEN
    IF DUP TYPE 3 ==                    @ There is object, check if array
    THEN
      DUP 'XY' STO SIZE 1 GET 'N' STO   @ Save Array and Number of points
      CLLCD                             @ Clear screen for progess updates
      LijCAL                            @ Calculate L coefficients
      KCAL                              @ Calculate final K coefficients
    ELSE EMESS                          @ Display error message
    END
    {CST POLY VIEWK Fx}                 @ Restore main variable order
    ORDER                               @ after new variables are created
  ELSE
    EMESS                               @ Display error message
  END
\>>
CST
{                                       @ Custom Menu for Polyfit program
  {"POLY" POLY}                         @ Main program to find coefficients
  {"VIEW" VIEWK}                        @ View polynomial coefficients
  {"Fx" Fx}                             @ Given X (on stack) return f(X)
  {}                                    @ Nothing
  {}                                    @ Nothing
  { "X-Y" XY}                           @ Recall input array to stack for editing
}
Cab
\<< @ Generate terms to sum for C(a)(b) @
  IF DUP2 ==
  THEN DROP2 1                          @ If a=b, the Cab=1
  ELSE
    0 1 \-> a b SUM Index
    \<<
      { 1 }                             @ Initial list

      DO  @ For all possible terms @
        WHILE @ Build list up to a-b terms @
          Index a b - <                 @ Until we hit last term
        REPEAT
          DUP Index GET 1 + 1 \->LIST +  @ Add next term to list
          Index 1 + 'Index' STO
        END

        DO  @ Increment last terms until limit @
          FETCH                         @ Get X values and multiply
          SUM + 'SUM' STO               @ Add to sum and save
          DUP Index GET 1 + Index SWAP
          DUP 4 ROLLD PUT               @ Increment last term and save
        UNTIL SWAP b Index + >          @ Until we go over b+Index limit
        END

        IF Index 1 >  @ Increment previous terms @
        THEN                            @ If possible
          DO
            Index 1 - 'Index' STO       @ Decrement index
            1 Index SUB                 @ Shorten list
            DUP Index GET 1 + Index SWAP
            DUP 4 ROLLD PUT             @ Increment this term and save
          UNTIL
            SWAP b Index + \<=          @ This term at limit
            Index 1 == OR               @ or this is last term
          END
        END
      UNTIL
        DUP Index GET b Index + >
        Index 1 == AND
      END
      DROP                              @ Remove last list from stack
      SUM
      -1 a b - ^ *                      @ Fix sign of product
    \>>
  END                                   @ If a=b
\>>
EMESS
\<< @ Display error message if no stack value, or value is not an array @
  CLLCD
  "Input array not found." 1 DISP
  "Use Matrix Writer to " 3 DISP
  "enter your X Y values" 4 DISP
  "into an array on the" 5 DISP
  "stack.  Then start" 6 DISP
  "POLY again." 7 DISP
  7 FREEZE
\>>
EQ
Y.EQ                                    @ For PLOTR to plot f(x)
FETCH
\<< @ Use Index list on stack to fetch X values and multiply @
  DUP SIZE 1 \-> lst n product          @ Save term list, size, and product
  \<<
    XY
    1 n
    FOR i
      DUP lst i GET 1 XYindex GET       @ Get Xi term
      product * 'product' STO           @ Multiply and save result
    NEXT
    DROP                                @ Drop XY array
    lst product                         @ Restore term list and product
  \>>
\>>
Fx
\<<  @ Calculate F(x) using K coefficients and x value from stack @
  K N GET                               @ First coefficient
  N 1 - 1
  FOR i                                 @ For i = N-1 to 1
    OVER *                              @ x times
    K i GET +                           @ Add in next coefficient
    -1
  STEP
  SWAP DROP                             @ Remove x value
\>>
KCAL
\<< @ Calculate K coefficients (f(x) = Kn*x^n-1 + ... K2*x + K1) @
  1 N
  START                                 @ Create K array with N zeroes
    0
  NEXT
  N \->ARRY 'K' STO                     @ Save K array

  1 N
  FOR i                                 @ i=1 to N
    "DOING K" STD i \->STR + 1 DISP     @ Display Ki value
    0                                   @ Zero sum
    i 1 - N 1 -
    FOR j                               @ j=i-1 to N-1
      @ Ki = Ki + Lj1 * C(j)(i-1) @
      L j 1 Lindex GET                  @ Get L term from array
      j i 1 - Cab                       @ Get Cab coefficient
      * +                               @ Sum to Ki
    NEXT
    K i 3 ROLL PUT 'K' STO              @ Add K value to array
    "K" i \->STR + "=" + K i GET \->STR + 3 DISP
  NEXT
\>>
LijCAL
@ Calculate triangular table of L values @
\<<
  "DOING L:" 1 DISP
  1 TMAX                                @ Create L values array
  START                                 @ With TMAX zeroes
    0
  NEXT
  TMAX \->ARRY                          @ L array initialized
  1 N                                   @ Set Lj,1 = Yj
  FOR j                                 @ For i=0, j = 1 to N
                                        @ L array already on stack
    j                                   @ L index value to place Y value
    XY j 2 * GET                        @ Y value for L array
    PUT                                 @ Y value now in L array
  NEXT
                                        @ Initialized L array left on stack
  1 N 1 -
  FOR i                                 @ For i = 1 to N-1
    1 N i -
    FOR j                               @ For j = 1 to N-i
      DUP DUP i 1 - j 1 + Lindex GET    @ Get L i-1,j+1 value
      SWAP i 1 - j Lindex GET -         @ Get L i-1,j value and subtract
      XY DUP i j + 1 XYindex GET        @ Get X i+j value
      SWAP j 1 XYindex GET -            @ Get X j value and subtract
      / i j Lindex SWAP PUT             @ Save new L i,j into array
    NEXT
  NEXT
  'L' STO                               @ Save final array
\>>
Lindex
@ Return index into L array of i,j element @
\<<
  \-> i j                               @ Save index values
  \<<
    i 1 - N i 2 / - * N + j +           @ index=(i-1)(N-i/2)+j+N
  \>>
\>>
PPAR
{
  (-6.5,-3.1)                           @ X range
  (6.5,3.2)                             @ Y range
  X                                     @ Independent Variable
  0                                     @ Resolution (pixel in each column)
  (0,0)                                 @ Axes intersection
  FUNCTION                              @ Function type
  Y                                     @ Dependent Variable
}
TMAX
@ Maximum Number of elements in the L triangular array @
\<<
  N DUP 1 + * 2 /                       @ TMAX = N(N+1)/2
\>>
VIEWK
\<< @ Display the K coefficients of the resulting polynomial expression @
  STD CLLCD "Hit VIEW to cont..." 1 DISP
  N 1
  FOR i                                 @ For each coefficient
    "K" i \->STR + "=" +
    K i GET \->STR + "*x^" + i 1 -
    \->STR + 3 DISP 0 WAIT DROP
    -1
  STEP
  CLLCD                                 @ Clear display to let user know
                                        @ the program is still working
  KILL                                  @ Get rid of HALT annunciator
\>>
XYindex
@ Return index into XY array of i,j element @
\<<
  \-> i j                               @ Save index values
  \<<
    i 1 - 2 * j +                       @ index=2(i-1)+j
  \>>
\>>
Y.EQ
\<< @ Fx function requires X value from stack, not compatible with PLOTR. @
  @ This program takes 'X' from PLOTR function, calls Fx, and @
  @ returns with f(x) value to satisfy PLOTR requirements. @
  X Fx
\>>
END_RPL